Gazelle distance sampling analysis for Tost 2025 data

knitr::opts_chunk$set(results=F, warning=F, message=F, fig.align='center', cache.lazy=F)

# Load some packages
library(plyr)
library(sf)
library(terra)
library(elevatr)
library(nngeo)
library(knitr)
library(ggplot2)
library(kableExtra)
library(leaflet)
library(RColorBrewer)
library(Distance)
library(glmmTMB)
library(patchwork)
library(rphylopic)
library(dplyr)
library(viridis)
library(readxl)

1. Introduction

In this notebook we analyse line transect data collected in the steppe habitats of Tost, Toson Bumba (Tost) Nature Reserve in September 2025. Detections were generated between 13-16 September 2025 along main 28 transects of variable length (ranging from 20 to 82 km) and undertaken by motorbike with single (14 transects) or double (14 transects) observer.

Here we apply a distance sampling approach to these data and generate abundance estimates for each survey. We will try to calculate the impact of environmental variables.

A one species were detected during the surveys - the goitered gazelle (Gazella subgutturosa)

The specific objectives of this notebook are therefore to:

• Estimate the detection function for goitered gazelles by single or double observer, NDVI, TRI, Water, HDI etc…
• To generate density and abundance estimates for the goitered gazelle in the survey 2025

2. Load and prepare the data

We begin by loading the csv files and spatial data.

# Set working directory
setwd("D:/SLCF/Ungulate/Steppe/Gazelle_DS_2025_Tost")
knitr::opts_knit$set(root.dir="D:/SLCF/Ungulate/Steppe/Gazelle_DS_2025_Tost")

# Survey data
Tost25 <- read.csv("D:/SLCF/Ungulate/Steppe/Gazelle_DS_2025_Tost/Tost_Sungulate_data2025.csv", check.names=F)

# Put them in a list chronologically
dat <- list(Tost25)

# Transects
#transectsG25 <- st_read("D:/SLCF/Ungulate/Steppe/Steppeshape/GGN_Sungulate_track2025.shp", quiet=T)
transectsT25 <- st_read("D:/SLCF/Ungulate/Steppe/Steppeshape/Tost_Sungulate_track2025.shp", quiet=T)

# Reproject transects to UTM 47N (best for Tost; 48N technically better for GGNP)
study.crs <- "EPSG:32647"

transectsT25 <- st_transform(transectsT25, crs=study.crs)
#transectsT25 <- st_transform(transectsT25, crs=study.crs)
study.crs <- "EPSG:32647"
# Define the survey area
# We can use an arbitrary 9km(GGN), 8km (Tost) buffer around MCP of transects, on the assumption that the habitat is the same within 10 km of transects
#ggn.area <- st_union(st_buffer(transectsG25, 9000))
tost.area <- st_union(st_buffer(transectsT25, 8000))
survey.area <- st_union(tost.area)
plot(survey.area)

# Note: as the survey area, we could also use the protected area boundaries (I think this is justifiable for GGN, We used MINUS/PLUS sampling).

# Also create a layer for the bounding box around this MCP
study.extent <- st_as_sfc(st_bbox(survey.area))

# Get elevation data for the survey area
elev <- get_elev_raster(st_as_sf(survey.area), z=9, clip="bbox") #110 m resolution
tri <- terrain(elev, v="TRI") # 3 x 3 window
#st_write(tost.area, "D:/SLCF/Ungulate/Steppe/Gazelle_DS_2025_Tost/TostFull_area.gpkg")

We then merge the data into one data-frame with standardised column names for easier analysis. In the process, we also make some minor fixes (removing blank rows, making the species names consistent etc.).

# Merge the data across sites and years

## i. Remove any blank rows
dat <- llply(dat, function(x) {
  x[x$Transect!="", ]
})
## ii. Remove some columns (which are blank anyway in the GGN data)
to.remove <- c("Leader", "", "Updated_by", "Category2", "IIIIII", "Updated_da", "Patrol_ID")
dat <- llply(dat, function(x) {
  x[, !names(x) %in% to.remove]
})
## iii. Add year and site columns
year.vec <- 2025
site.vec <- c("Tost")
for (i in 1:length(dat)) {
dat[[i]]$Year1 <- year.vec[i]
dat[[i]]$Site <- site.vec[i]}

## iv. Remove all underscores and spaces from column names
dat <- llply(dat, function(x) {
  names(x) <- tolower(names(x))
  #names(x) <- gsub("\\.", "", names(x)) #for removing periods
  names(x) <- gsub("_", "", names(x)) 
  names(x) <- gsub(" ", "", names(x)) 
  x
})


## v. Collapse into one data-frame
dat <- ldply(dat, function(x) {x}, .id=NULL)

# vi. Remove all rows with no species information. This includes rows for transects with no detections (these transects are captured in the spatial data anyway). 
dat <- dat[dat$species!="", ]

# vii. Standardise the species names
dat$species <- factor(dat$species)
levels(dat$species) <- c( "goitered_gazelle")

Next we process the spatial data. First, we clean up and merge all the transect data.

# Rename the transect column to just 'name'
#names(transectsG25)[names(transectsG25)=="Track_Firs"] <- "name"
names(transectsT25)[names(transectsT25)=="Track_Firs"] <- "name"

# Create a standardised column called "matchID", which will be a combination of the transect name and year
#transectsG25$matchID <- paste0(transectsG25$name, "-2025")
transectsT25$matchID <- paste0(transectsT25$name, "-2025")

# Also create a "surveyID" column, i.e. a combination of the site name and year
#transectsG25$surveyID <- paste0("GGN", "-2025")
transectsT25$surveyID <- paste0("Tost", "-2025")

# Combine the transects into one layer
transects <- rbind(transectsT25[, c("name", "matchID", "surveyID")])

# Add the length of each transect (in km)
transects$length <- as.numeric(st_length(transects) / 1000)
# Create the same matchID and surveyID columns in the detection data
dat$matchID <- paste0(dat$transect, "-", dat$year)
dat$surveyID <- paste0(dat$site, "-", dat$year)

# Lump together motorbike
dat$vehicle <- factor(dat$transectm)
levels(dat$vehicle) <- c("moto")
# Lump together Single and Double
dat$observern <- factor(dat$observern)
levels(dat$observern) <- c("double","single")
# Lump together 0 and 1
dat$foresttt <- factor(dat$foresttt)
levels(dat$foresttt) <- c("0", "1")
# Lump together 0 and 1
dat$watertt <- factor(dat$watertt)
levels(dat$watertt) <- c("0", "1")
# Lump together ndvi categories
dat$ndvitt <- factor(dat$ndvitt)
levels(dat$ndvitt) <- c("0.05-0.10", "0.11-0.15", "0.16-0.20")
# Lump together hdi categories
dat$hditt <- factor(dat$hditt)
levels(dat$hditt) <- c("0.00-0.01", "0.02-0.03", "0.04-0.05", ">0.06")
# Lump together tri categories
dat$tritt <- factor(dat$tritt)
levels(dat$tritt) <- c("1.41-2.00", "3.00-4.00", "5.00-6.00", ">7")
# Check that all the rows in the detection data have tracks associated with them
#all(dat$matchID %in% transects$matchID) #TRUE

Next we calculate the perpendicular distance of each detection from the transect line. These distances are the critical data input to distance sampling.

# We first simplify the transect tracks a bit, because we want the perpendicular distance from the main direction of travel of the motorbike, discounting the minor turns and wiggles in the track. 
# We use a threshold distance of 100 m for the simplification of the GPS tracks, which preserves the overall structure.

transects.simple <- st_simplify(transects, dTolerance = 100)
# Reproject the lat longs to UTM 47N
observ.locations <- st_as_sf(dat, coords=c("longitude", "latitude"), crs="EPSG:4326")
observ.locations <- st_transform(observ.locations, crs=study.crs)
dat <- cbind(dat, as.data.frame(st_coordinates(observ.locations)))

# Calculate the object XY coordinates

dat$objectX <- dat$X + (sin(dat$angle * (pi/180)) * dat$distance)
dat$objectY <- dat$Y + (cos(dat$angle * (pi/180)) * dat$distance)

# We could try and use the near distance and near lat-longs, but still we don't know which side of the transect the animals were on, so we still can't pinpoint their coordinates. For these, we will just use the near distance measurements (presumably measured in the field).

# Create a layer for all the detections with angles/distances
object.locations <- st_as_sf(dat[!is.na(dat$objectX), ], coords=c("objectX", "objectY"), crs=study.crs)

# Calculate the perpendicular distances to the simplified transect lines
object.locations <- adply(object.locations, 1, function(x){
  transect.x <- transects.simple[transects.simple$matchID==x$matchID, ]
  dist.x <- st_nn(x, transect.x, returnDist=T, progress=F)
  unlist(dist.x$dist)
})
names(object.locations)[names(object.locations)=="V1"] <- "perpdist"

plot(object.locations$perpdist, object.locations$perpdist) #note; not need

#plot(object.locations$minimumdistance, object.locations$perpdist)
#Plot the relationship between the field-estimated distances and the perpendicular distances calculated from the angle/distance data
#plot(object.locations$minimumdistance, object.locations$perpdist)
# There is a positive relationship, but a lot of scatter. Why would there be quite big differences in some cases?

# Also output the connecting lines for visualising in QGIS
perpdist.lines <- adply(object.locations, 1, function(x){
  transect.x <- transects.simple[transects.simple$name==x$transect, ]
  dist.x <- st_nn(x, transect.x, progress=F)
  line.x <- st_connect(x, transect.x, ids=dist.x, dist=10)
  st_as_sf(line.x)
})
st_geometry(perpdist.lines) <- perpdist.lines$x

#st_write(observ.locations,"D:/SLCF/Ungulate/Steppe/Gazelle_DS_2025_Tost/Tost_observer_locations_2025.gpkg")
#st_write(object.locations, "D:/SLCF/Ungulate/Steppe/Gazelle_DS_2025_Tost/Tost_object_locations_2025.gpkg")
#st_write(transects.simple, "D:/SLCF/Ungulate/Steppe/Gazelle_DS_2025_Tost/Tost_tracks_2025_simple_100.gpkg")
#st_write(perpdist.lines, "D:/SLCF/Ungulate/Steppe/Gazelle_DS_2025_Tost/Tost_perpendicular_distances_2025.gpkg")

# Now add the perpendicular distances into the detection data
dat$perpdist <- NA
dat[!is.na(dat$objectX), "perpdist"] <- as.data.frame(object.locations)[, "perpdist"]

# Fill in the perpendicular distances for the rows with missing distance/angle data using the near distances recorded in the field
dat$perpdist[is.na(dat$perpdist)] <- dat$minimumdistance[is.na(dat$perpdist)]

#write.csv(dat, file="Tost_steppe_surveydata_2025_cleaned_V4.csv", row.names=F)
#write.csv(dat, file="steppe_surveydata_2025_cleaned.csv", row.names=F)
# Create a data subset using the goitered gazelle detections only
dat.ggaz <- dat[dat$species=="goitered_gazelle", ]

As the final preparation step, we create the three data-frames required by Distance, which are: i) a region table, identifying any spatial strata), ii) a sample table, providing information on the transects, and iii) an observation table, providing the distance data.

Note: it’s also possible to provide Distance with one table with all of these data (in a so-called ‘flatfile’), but this is a less flexible approach. For example, we can’t customise the survey areas using this approach.

# Create the standard objects (region.table, sample.table and obs.table) required for conventional distance sampling. 

# The region table needs columns 'Region.Label' and 'Area' (if abundance estimates are required, otherwise it can be set to 0)

#region.table <- data.frame(Region.Label=unique(dat.ggaz$surveyID), Area=as.numeric(c(st_area(tost.area), rep(st_area(ggn.area), 2)) / 1000^2))
region.table <- data.frame(Region.Label=unique(dat.ggaz$surveyID), Area=as.numeric(c(st_area(tost.area)) / 1000^2))

# The sample table needs columns 'Sample.Label', 'Region.Label' and 'Effort'
sample.table <- data.frame(Sample.Label=transects.simple$matchID, Region.Label=transects.simple$surveyID, Effort=transects$length) #note, transect length is in kilometres

obs.table <- data.frame(object=1:nrow(dat.ggaz), distance=dat.ggaz$perpdist, size=dat.ggaz$groupsize, Region.Label=dat.ggaz$surveyID, Sample.Label=dat.ggaz$matchID, observern=dat.ggaz$observern, hditt=dat.ggaz$hditt, foresttt=dat.ggaz$foresttt, tritt=dat.ggaz$tritt, watertt=dat.ggaz$watertt, ndvitt=dat.ggaz$ndvitt)

3. Initial exploration of the data

We’ll now explore the contents of the data a bit.

First let’s calculate sample sizes for each species in each survey.

3.1. Sample sizes

kable(as.data.frame.array(table(dat$species, factor(dat$year1))), col.names=c("Species", "Tost 2025")) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width=F)

Species	Tost 2025
goitered_gazelle	93

We can see that there are few detections of khulan. The detections are overwhelmingly of goitered gazelle.

Next let’s create maps of the transects and animal detections for the survey.

3.2. Mapping the detections

species.pal <- colorFactor(palette=brewer.pal(3, "Set1"), domain=object.locations$species)
leaflet() %>% 
  addTiles(urlTemplate = "http://mt1.google.com/vt/lyrs=s&x={x}&y={y}&z={z}") %>% 
  addPolylines(data=st_transform(survey.area, "EPSG:4326"), color="red", weight=1.8, opacity=0.6, dashArray="5") %>%
  addPolylines(data=st_transform(st_zm(transects.simple, drop=T, what="ZM"), "EPSG:4326"), color="black", weight=1.4, opacity=0.5) %>%
  addPolylines(data=st_transform(perpdist.lines, "EPSG:4326"), color="white", weight=1.3, opacity=0.5) %>%
  addCircleMarkers(data=st_transform(object.locations, "EPSG:4326"), radius=4, stroke=F, fillOpacity=0.9, color=~species.pal(object.locations$species)) %>% 
  addLegend("topright", pal=species.pal, values=object.locations$species,
            title = "Species location")

Goitered gazelle appear to be widely distributed across the survey area, except in the plains to west of the Tost, where few detections were made.

Next we plot the detection distances for the species in survey.

3.3. Distance histograms

# Plot a histogram of the distances
ggplot(dat, aes(x=perpdist, y=after_stat(density))) +
  geom_histogram(bins=30, fill="grey", alpha=0.8) +
  geom_density(colour="darkred", alpha=0.8) +
  labs(x="Distance from transect (m)", subtitle="Detection distances by species and site") +
  theme_minimal() +
  facet_wrap(~species+site+year, ncol=2) 

The goitered gazelle data show the expected decline in detection rate with distance. However, in Tost (2025), there seems to be a high of detection at very close distance. Farther more 250m very lack detection. If this is the case, density estimates will be negatively biased. To attempt to account for this, we could bin the data in the first 100-150 m, in effect making the assumption that animals pushed out from 0 m are nonetheless caught somewhere in the range 0-150 m. The other option would be to apply 2D distance sampling (we could explore this in a separate notebook).

For gazelle data, there seems to be a slight excess of detections at ~750 m. Standard practice is to truncate the furthest 5-10% of the distances, which would mean truncating at 746 (90% quantile) or 907 (95% quantile). Given the potential outliers at ~750 m, we will truncate at 10 %.

3.4. Data subsetted by observer type

We will now calculate sample sizes for the two observer types used in the survey. We should also inspect histograms of the distances separately for both observer type.

We’ll do this for the goitered gazelle data only, since these are the data we will be focusing on.

kable(as.data.frame.array(table(dat.ggaz$observern)), col.names=c("Observer type", "n")) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width=F)

Observer type	n
double	53
single	40

There are perhaps sufficient detection to fit separate detection functions for double and single. Its 14=14 equal in Tost. For Tost, 40 sightings were observed by single observer and 53 sightings observed by double observer type. I think it is not too much difference.

# Plot a histogram of the distances
ggplot(dat.ggaz, aes(x=perpdist, y=after_stat(density))) +
  geom_histogram(bins=30, fill="grey", alpha=0.8) +
  geom_density(colour="darkred", alpha=0.8) +
  labs(x="Distance from transect (m)", subtitle="Goitered gazelle distances by observer type") +
  theme_minimal() +
  facet_wrap(~observern, ncol=2) 

In the different observer type, the highest density of observations is made in the first 200m, but the single observer’s observations are higher in close range compared to double observer. However, the double observer shows higher observations density at distances beyond 350m compared to the single observer.

# Plot a histogram of the distances
ggplot(dat.ggaz, aes(x=perpdist, y=after_stat(density))) +
  geom_histogram(bins=30, fill="grey", alpha=0.8) +
  geom_density(colour="darkred", alpha=0.8) +
  labs(x="Distance from transect (m)", subtitle="Goitered gazelle distances by shrub or forest(0-absence 1-presence)") +
  theme_minimal() +
  facet_wrap(~foresttt, ncol=2) 

kable(as.data.frame.array(table(dat.ggaz$foresttt)), col.names=c("Shrub 0-absence 1-presence", "n")) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width=F)

It appears that the influence of the forest has a significant impact on detection in terms of distance and density.

# Plot a histogram of the distances
ggplot(dat.ggaz, aes(x=perpdist, y=after_stat(density))) +
  geom_histogram(bins=30, fill="grey", alpha=0.8) +
  geom_density(colour="darkred", alpha=0.8) +
  labs(x="Distance from transect (m)", subtitle="Goitered gazelle distances by water (0-absence 1-presence)") +
  theme_minimal() +
  facet_wrap(~watertt, ncol=2) 

kable(as.data.frame.array(table(dat.ggaz$watertt)), col.names=c("Water 0-absence 1-presence", "n")) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width=F)

# Plot a histogram of the distances
ggplot(dat.ggaz, aes(x=perpdist, y=after_stat(density))) +
  geom_histogram(bins=30, fill="grey", alpha=0.8) +
  geom_density(colour="darkred", alpha=0.8) +
  labs(x="Distance from transect (m)", subtitle="Goitered gazelle distances by TRI(categories)") +
  theme_minimal() +
  facet_wrap(~tritt, ncol=2) 

kable(as.data.frame.array(table(dat.ggaz$tritt)), col.names=c("TRI categories", "n")) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width=F)

# Plot a histogram of the distances
ggplot(dat.ggaz, aes(x=perpdist, y=after_stat(density))) +
  geom_histogram(bins=30, fill="grey", alpha=0.8) +
  geom_density(colour="darkred", alpha=0.8) +
  labs(x="Distance from transect (m)", subtitle="Goitered gazelle distances by NDVI (categories)") +
  theme_minimal() +
  facet_wrap(~ndvitt, ncol=2) 

kable(as.data.frame.array(table(dat.ggaz$ndvitt)), col.names=c("NDVI categories", "n")) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width=F)

# Plot a histogram of the distances
ggplot(dat.ggaz, aes(x=perpdist, y=after_stat(density))) +
  geom_histogram(bins=30, fill="grey", alpha=0.8) +
  geom_density(colour="darkred", alpha=0.8) +
  labs(x="Distance from transect (m)", subtitle="Goitered gazelle distances by HDI (categories)") +
  theme_minimal() +
  facet_wrap(~hditt, ncol=2) 

kable(as.data.frame.array(table(dat.ggaz$hditt)), col.names=c("HDI categories", "n")) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width=F)

For gazelles, proximity/near (5000m) to a water source seems to be of little importance.

3.5. Size-biased detection

To plot herd size as a function of distance, testing the hypothesis that herds detected at far distances are larger than those at close distances. This is important because we might be missing smaller herds at further distances, and we therefore risk producing negatively biased abundance estimates if we do not correct for this effect.

Note, we will also introduce herd size as a potential covariate when we are modelling the detection function to test for this effect. But it is nonetheless useful also to plot the data as we do here.

ggplot(dat.ggaz, aes(y=groupsize, x=perpdist)) +
  geom_point() +
  geom_smooth(se=T) +
  labs(x="Distance (m)", y="Gazelle herd size") +
  theme_minimal()

There is perhaps a bit tendency for herd size to increase at further distances, implying we might be missing some small groups, but this relationship doesn’t look to be statistically clear. But there is not much difference within 1000 meters.

We can fit a zero-truncated Poisson model to investigate if the relationship between herd size and distance is statistically clear. The Poisson distribution is appropriate because herd sizes are integers and cannot be negative. The zero-truncation is necessary because herd sizes of zero cannot be observed (in practice, a Poisson generalized linear model will also give very similar results). We will truncate the furthest 10% of the distances.

This figure shows that the human impact index is highly dependent on herd size and detection.

# Truncation distance
trunc.dist <- round(quantile(dat.ggaz$perpdist, 0.9), 0)


# Zero-truncated Poisson model of herd size as a function of distance
mod.size <- glmmTMB(groupsize ~ perpdist, family=list(family="truncated_poisson", link="log"), data=dat.ggaz[dat.ggaz$perpdist < trunc.dist, ])
summary(mod.size)

##  Family: truncated_poisson  ( log )
## Formula:          groupsize ~ perpdist
## Data: dat.ggaz[dat.ggaz$perpdist < trunc.dist, ]
## 
##      AIC      BIC   logLik deviance df.resid 
##    726.2    731.0   -361.1    722.2       81 
## 
## 
## Conditional model:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 1.8616989  0.0613668  30.337   <2e-16 ***
## perpdist    0.0002120  0.0002223   0.954     0.34    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat("\nConfidence intervals:\n")

## 
## Confidence intervals:

confint(mod.size)

##                     2.5 %       97.5 %     Estimate
## (Intercept)  1.7414221665 1.9819756994 1.8616989329
## perpdist    -0.0002236292 0.0006476505 0.0002120106

# For comparison, we can also fit a Poisson generalized linear model
#mod.size.glm <- glm(groupsize ~ perpdist, family="poisson", data=dat.ggaz[dat.ggaz$perpdist < trunc.dist, ])
#summary(mod.size.glm) # slope estimate is very simlar to the zero-truncated model

We can see that there does seem to be a size-bias, with the slope being significant different from zero (p < 0.05). The confidence intervals for the slope estimate also do not cross zero.

We can plot the model predictions.

# Make predictions of group size according to distance
mod.size.pred <- predict(mod.size, newdata=data.frame(perpdist=0:trunc.dist), type="link", se=T)

# Put the data together for plotting
mod.size.pred <- data.frame(Distance=0:trunc.dist, groupsize=exp(mod.size.pred$fit), LCI=exp(mod.size.pred$fit - 1.96 * mod.size.pred$se.fit), UCI=exp(mod.size.pred$fit + 1.96 * mod.size.pred$se.fit))

# Plot it
ggplot(data=mod.size.pred, aes(x=Distance, y=groupsize)) +
  geom_line() +
  geom_ribbon(aes(ymin=LCI, ymax=UCI), alpha=0.25) +
  lims(y=c(0, 20)) +
  labs(y="Herd size", x="Distance (m)") +
  theme_minimal()

The average herd size in the observed data was 6.71, whilst our model predicts that group sizes were actually 6.43 (using the predicted group size at 0 m, i.e. assuming we perfectly observe group size at zero distance).

4. Identifying an appropriate detection function

We will now formally model the detection process using Distance. For now, we will not look at the abundance estimates, and instead just focus on identifying a detection function which parsimoniously explains the data.

Our methodology will be as follows:

1. Define the most complicated model (‘maximal model’) and test goodness-of-fit

2. Test if goodness-of-fit improves if the data < 150 m are binned

3. Carry out model selection using AIC

4. Select the top model (lowest AIC) for making predictions

Comparing the goodness-of-fit in (1) and (2) will help to identify if animals are being pushed out from the transect line and if this is affecting the results.

For model selection, we will fit null models (with no covariates) and also models with observer type as a covariate. This will allow the scale parameter of the detection function to vary by observer type. We will also fit herd size, to account for the possibility of size-biased detection.

Regarding the type of detection function, we will try both the half-normal and hazard detection functions (both commonly used in distance sampling), with and without adjustment terms. We will let Distance automatically find the appropriate number of adjustment terms.

4.1. Goodness-of-fit of continuous and binned data

dsmod.hr.max.obs <- ds(obs.table, truncation="10%", formula=~observern+size,  key="hr", adjustment=NULL, cutpoints=NULL, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)

(gof.max <- gof_ds(dsmod.hr.max.obs))

## 
## Goodness of fit results for ddf object
## 
## Distance sampling Cramer-von Mises test (unweighted)
## Test statistic = 0.0420581 p-value = 0.922259

#Distance sampling Cramer-von Mises test (unweighted)
#Test statistic =

The model fits the data OK (p > 0.05). However, it does seem like there are ‘too few’ observations at close distances. The problem does not look too serious though. ##Goodness of fit results for ddf object

##Distance sampling Cramer-von Mises test (unweighted) ##Test statistic = 0.0420581 p-value = 0.922259

# The most complicated model will include covariates for observer type and herd size.
# The hazard-rate is more complicated than the half-normal (2 parameters instead of 1), so we will use it.
dsmod.hr.max.hdi <- ds(obs.table, truncation="10%", formula=~hditt+size,  key="hr", adjustment=NULL, cutpoints=NULL, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)

(gof.max <- gof_ds(dsmod.hr.max.hdi))

#Distance sampling Cramer-von Mises test (unweighted)

The model fits the data OK (p > 0.05). However, it does seem like there are ‘too few’ observations at close distances. The problem does not look too serious though. ##Goodness of fit results for ddf object

##Distance sampling Cramer-von Mises test (unweighted) ##Test statistic = 0.0256372 p-value = 0.988217

dsmod.hr.max.forest <- ds(obs.table, truncation="10%", formula=~foresttt+size,  key="hr", adjustment=NULL, cutpoints=NULL, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)

(gof.max <- gof_ds(dsmod.hr.max.forest))

#Distance sampling Cramer-von Mises test (unweighted)
#Test statistic = 

dsmod.hr.max.tri <- ds(obs.table, truncation="10%", formula=~tritt+size,  key="hr", adjustment=NULL, cutpoints=NULL, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)

(gof.max <- gof_ds(dsmod.hr.max.tri))

#Distance sampling Cramer-von Mises test (unweighted)
#Test statistic =

dsmod.hr.max.water <- ds(obs.table, truncation="10%", formula=~watertt+size,  key="hr", adjustment=NULL, cutpoints=NULL, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)

(gof.max <- gof_ds(dsmod.hr.max.water))

#Distance sampling Cramer-von Mises test (unweighted)
#Test statistic = 

dsmod.hr.max.ndvi <- ds(obs.table, truncation="10%", formula=~ndvitt+size,  key="hr", adjustment=NULL, cutpoints=NULL, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)

(gof.max <- gof_ds(dsmod.hr.max.ndvi))

#Distance sampling Cramer-von Mises test (unweighted)
#Test statistic = 

Goodness of fit results for ddf object #Distance sampling Cramer-von Mises test (unweighted) #Test statistic = 0.0420581 p-value = 0.922259 obs

#Distance sampling Cramer-von Mises test (unweighted) #Test statistic = 0.0259809 p-value = 0.987436 hdi

#Distance sampling Cramer-von Mises test (unweighted) #Test statistic = 0.014209 p-value = 0.999763 forest

#Distance sampling Cramer-von Mises test (unweighted) #Test statistic = 0.0346663 p-value = 0.958446 tri

#Distance sampling Cramer-von Mises test (unweighted) #Test statistic = 0.0242996 p-value = 0.990977 water

#Distance sampling Cramer-von Mises test (unweighted) #Test statistic = 0.0235485 p-value = 0.992336 ndvi

Let’s test goodness-of-fit with binned data. Goodness of fit results for ddf object

Chi-square tests

#P = 0.84836 with 18 degrees of freedom observer #P = 0.67726 with 16 degrees of freedom hdi #P = 0.81551 with 18 degrees of freedom forest #P = 0.67373 with 16 degrees of freedom TRI #P = 0.81278 with 18 degrees of freedom water #P = 0.76043 with 17 degrees of freedom NDVI

ds.cutpoints <- c(0, 50, seq(75, trunc.dist, 25))

dsmod.hr.max.obs.binned <- ds(obs.table, formula=~observern+size, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)

(gof.max.obs.binned <- gof_ds(dsmod.hr.max.obs.binned))

ds.cutpoints <- c(0, 50, seq(75, trunc.dist, 25))

dsmod.hr.max.hdi.binned <- ds(obs.table, formula=~hditt+size, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)

(gof.max.hdi.binned <- gof_ds(dsmod.hr.max.hdi.binned))

ds.cutpoints <- c(0, 50, seq(75, trunc.dist, 25))

dsmod.hr.max.forest.binned <- ds(obs.table, formula=~foresttt+size, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)

(gof.max.forest.binned <- gof_ds(dsmod.hr.max.forest.binned))

ds.cutpoints <- c(0, 50, seq(75, trunc.dist, 25))

dsmod.hr.max.tri.binned <- ds(obs.table, formula=~tritt+size, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)

(gof.max.tri.binned <- gof_ds(dsmod.hr.max.tri.binned))

ds.cutpoints <- c(0, 50, seq(75, trunc.dist, 25))

dsmod.hr.max.water.binned <- ds(obs.table, formula=~watertt+size, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)

(gof.max.water.binned <- gof_ds(dsmod.hr.max.water.binned))

ds.cutpoints <- c(0, 50, seq(75, trunc.dist, 25))

dsmod.hr.max.ndvi.binned <- ds(obs.table, formula=~ndvitt+size, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)

(gof.max.ndvi.binned <- gof_ds(dsmod.hr.max.ndvi.binned))

The binned (all)data also appear to fit the model OK. Note we don’t get a quantile-quantile plot when using binned data (it is only for continuous data). The number of gazelle expected in the covered region is 300.9. For the binned data, it is 286.95. The fact that the estimate for the binned data is slightly higher is in line with the expectation that we might be underestimating abundance due to the disturbance effect, but the effect is minimal.

We can plot the detection functions for the continuous and binned data to see how the detection functions compare.

plot(dsmod.hr.max.obs)

plot(dsmod.hr.max.obs.binned)

plot(dsmod.hr.max.hdi)

plot(dsmod.hr.max.hdi.binned)

plot(dsmod.hr.max.forest)

plot(dsmod.hr.max.forest.binned)

plot(dsmod.hr.max.tri)

plot(dsmod.hr.max.tri.binned)

plot(dsmod.hr.max.water)

plot(dsmod.hr.max.water.binned)

plot(dsmod.hr.max.ndvi)

plot(dsmod.hr.max.ndvi.binned)

Note: the points are the predicted detection probabilities for each observation. They do not lie exactly on the overall detection function line because the covariate type and herd sizes differ between observations. we will use the continuous data going forward.

4.2. Model selection

Next we will run a set of models and rank them according to AIC.

# Null models

## Half-norm
dsmod.hn.0 <- ds(obs.table, formula=~1, key="hn", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="strict", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)

## Hazard
dsmod.hr.0 <- ds(obs.table, formula=~1, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="strict", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)

# With adjustment terms (Distance automatically selects the best number of terms to use based on AIC)
dsmod.hn.adj <- ds(obs.table, formula=~1, key="hn", adjustment="cos", cutpoints=ds.cutpoints, monotonicity="strict", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
dsmod.hr.adj <- ds(obs.table, formula=~1, key="hr", adjustment="cos", cutpoints=ds.cutpoints, monotonicity="strict", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)

# Covariates

## HDI
dsmod.hn.1 <- ds(obs.table, formula=~hditt, key="hn", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
dsmod.hr.1 <- ds(obs.table, formula=~hditt, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)

## Herd size
dsmod.hn.2 <- ds(obs.table, formula=~size, key="hn", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
dsmod.hr.2 <- ds(obs.table, formula=~size, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
## observer type
dsmod.hn.3 <- ds(obs.table, formula=~observern, key="hn", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
dsmod.hr.3 <- ds(obs.table, formula=~observern, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
## Forest
dsmod.hn.4 <- ds(obs.table, formula=~foresttt, key="hn", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
dsmod.hr.4 <- ds(obs.table, formula=~foresttt, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)

# HDI-size
dsmod.hn.5 <- ds(obs.table, formula=~hditt+size, key="hn", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
dsmod.hr.5 <- ds(obs.table, formula=~hditt+size, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)

# Observer-size
dsmod.hn.6 <- ds(obs.table, formula=~observern+size, key="hn", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
dsmod.hr.6 <- ds(obs.table, formula=~observern+size, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)

# forest-size
dsmod.hn.7 <- ds(obs.table, formula=~foresttt+size, key="hn", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
dsmod.hr.7 <- ds(obs.table, formula=~foresttt+size, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
# Orbsever-size-hdi
dsmod.hn.8 <- ds(obs.table, formula=~observern+hditt+size, key="hn", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
dsmod.hr.8 <- ds(obs.table, formula=~observern+hditt+size, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
# Orbsever-size-forest
dsmod.hn.9 <- ds(obs.table, formula=~observern+foresttt+size, key="hn", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
dsmod.hr.9 <- ds(obs.table, formula=~observern+foresttt+size, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
# Orbsever-hdi-forest
dsmod.hn.10 <- ds(obs.table, formula=~observern+hditt+foresttt, key="hn", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
dsmod.hr.10 <- ds(obs.table, formula=~observern+hditt+foresttt, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
# hdi-forest
dsmod.hn.11 <- ds(obs.table, formula=~hditt+foresttt, key="hn", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
dsmod.hr.11 <- ds(obs.table, formula=~hditt+foresttt, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
# TRI
dsmod.hn.12 <- ds(obs.table, formula=~tritt, key="hn", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
dsmod.hr.12 <- ds(obs.table, formula=~tritt, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
# TRI-size
dsmod.hn.13 <- ds(obs.table, formula=~tritt+size, key="hn", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
dsmod.hr.13 <- ds(obs.table, formula=~tritt+size, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
# Observer-TRI-size
dsmod.hn.14 <- ds(obs.table, formula=~observern+tritt+size, key="hn", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
dsmod.hr.14 <- ds(obs.table, formula=~observern+tritt+size, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
# TRI-size-forest
dsmod.hn.15 <- ds(obs.table, formula=~tritt+foresttt+size, key="hn", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
dsmod.hr.15 <- ds(obs.table, formula=~tritt+foresttt+size, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
# Observer-TRI-size-forest
dsmod.hn.16 <- ds(obs.table, formula=~observern+tritt+foresttt+size, key="hn", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
dsmod.hr.16 <- ds(obs.table, formula=~observern+tritt+foresttt+size, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
# Observer-TRI-forest
dsmod.hn.17 <- ds(obs.table, formula=~observern+tritt+foresttt, key="hn", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
dsmod.hr.17 <- ds(obs.table, formula=~observern+tritt+foresttt, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
# Observer-TRI
dsmod.hn.18 <- ds(obs.table, formula=~observern+tritt, key="hn", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
dsmod.hr.18 <- ds(obs.table, formula=~observern+tritt, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
# Observer-TRI-size-forest-hdi
dsmod.hn.19 <- ds(obs.table, formula=~observern+tritt+foresttt+hditt+size, key="hn", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
dsmod.hr.19 <- ds(obs.table, formula=~observern+tritt+foresttt+hditt+size, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
# TRI-size-hdi
dsmod.hn.20 <- ds(obs.table, formula=~observern+tritt+hditt+size, key="hn", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
dsmod.hr.20 <- ds(obs.table, formula=~observern+tritt+hditt+size, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
# Water
dsmod.hn.21 <- ds(obs.table, formula=~watertt, key="hn", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
dsmod.hr.21 <- ds(obs.table, formula=~watertt, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
# Water-size
dsmod.hn.22 <- ds(obs.table, formula=~watertt+size, key="hn", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
dsmod.hr.22 <- ds(obs.table, formula=~watertt+size, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
# Water-tri-size
dsmod.hn.23 <- ds(obs.table, formula=~watertt+tritt+size, key="hn", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
dsmod.hr.23 <- ds(obs.table, formula=~watertt+tritt+size, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
# Water-forest-size
dsmod.hn.24 <- ds(obs.table, formula=~watertt+foresttt+size, key="hn", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
dsmod.hr.24 <- ds(obs.table, formula=~watertt+foresttt+size, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
# ndvi
dsmod.hn.25 <- ds(obs.table, formula=~ndvitt, key="hn", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
dsmod.hr.25 <- ds(obs.table, formula=~ndvitt, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
# ndvi-size
dsmod.hn.26 <- ds(obs.table, formula=~ndvitt+size, key="hn", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
dsmod.hr.26 <- ds(obs.table, formula=~ndvitt+size, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
# ndvi-tri-size
dsmod.hn.27 <- ds(obs.table, formula=~ndvitt+tritt+size, key="hn", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
dsmod.hr.27 <- ds(obs.table, formula=~ndvitt+tritt+size, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
# ndvi-forest-size
dsmod.hn.28 <- ds(obs.table, formula=~ndvitt+foresttt+size, key="hn", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
dsmod.hr.28 <- ds(obs.table, formula=~ndvitt+foresttt+size, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
# ndvi-water-size
dsmod.hn.29 <- ds(obs.table, formula=~ndvitt+watertt+size, key="hn", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)
dsmod.hr.29 <- ds(obs.table, formula=~ndvitt+watertt+size, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)

# Output a summary table

kbl(summarize_ds_models(dsmod.hn.0, dsmod.hr.0, dsmod.hn.adj, dsmod.hr.adj, dsmod.hn.1, dsmod.hr.1, dsmod.hn.2, dsmod.hr.2, dsmod.hn.3, dsmod.hr.3, dsmod.hn.4, dsmod.hr.4, dsmod.hn.5, dsmod.hr.5, dsmod.hn.6, dsmod.hr.6, dsmod.hn.7, dsmod.hr.7, dsmod.hn.8, dsmod.hr.8, dsmod.hn.9, dsmod.hr.9, dsmod.hn.10, dsmod.hr.10, dsmod.hn.11, dsmod.hr.11, dsmod.hn.12, dsmod.hr.12, dsmod.hn.13, dsmod.hr.13, dsmod.hn.14, dsmod.hr.14, dsmod.hn.15, dsmod.hr.15, dsmod.hn.16, dsmod.hr.16, dsmod.hn.17, dsmod.hr.17, dsmod.hn.18, dsmod.hr.18, dsmod.hn.19, dsmod.hr.19, dsmod.hn.20, dsmod.hr.20, dsmod.hn.21, dsmod.hr.21, dsmod.hn.22, dsmod.hr.22, dsmod.hn.23, dsmod.hr.23, dsmod.hn.24, dsmod.hr.24, dsmod.hn.25, dsmod.hr.25, dsmod.hn.26, dsmod.hr.26, dsmod.hn.27, dsmod.hr.27, dsmod.hn.28, dsmod.hr.28, dsmod.hr.29, dsmod.hn.29), digits=3) %>%
  kable_material(c("striped", "hover"), full_width=F) 

	Model	Key function	Formula	\(\chi^2\) \(p\)-value	\(\hat{P_a}\)	se(\(\hat{P_a}\))	\(\Delta\)AIC
24		Hazard-rate	~observern + hditt + foresttt	0.753	0.291	0.058	0.000
10		Hazard-rate	~observern	0.874	0.286	0.058	1.747
2		Hazard-rate	~1	0.918	0.279	0.057	1.859
4		Hazard-rate	~1	0.918	0.279	0.057	1.859
20		Hazard-rate	~observern + hditt + size	0.673	0.283	0.056	2.205
42		Hazard-rate	~observern + tritt + foresttt + hditt + size	0.504	0.291	0.045	2.773
46		Hazard-rate	~watertt	0.887	0.283	0.057	3.104
3		Half-normal with cosine adjustment term of order 2	~1	0.923	0.366	0.039	3.173
12		Hazard-rate	~foresttt	0.892	0.272	0.058	3.294
16		Hazard-rate	~observern + size	0.831	0.289	0.058	3.326
6		Hazard-rate	~hditt	0.773	0.291	0.056	3.553
22		Hazard-rate	~observern + foresttt + size	0.798	0.274	0.061	3.580
8		Hazard-rate	~size	0.891	0.280	0.057	3.827
54		Hazard-rate	~ndvitt	0.857	0.263	0.061	3.952
48		Hazard-rate	~watertt + size	0.854	0.284	0.057	5.097
14		Hazard-rate	~hditt + size	0.736	0.282	0.057	5.144
18		Hazard-rate	~foresttt + size	0.861	0.268	0.059	5.207
26		Hazard-rate	~hditt + foresttt	0.727	0.289	0.057	5.380
56		Hazard-rate	~ndvitt + size	0.817	0.264	0.061	5.891
44		Hazard-rate	~observern + tritt + hditt + size	0.458	0.287	0.051	6.216
52		Hazard-rate	~watertt + foresttt + size	0.820	0.272	0.058	6.556
38		Hazard-rate	~observern + tritt + foresttt	0.686	0.280	0.053	6.661
40		Hazard-rate	~observern + tritt	0.731	0.292	0.052	7.261
61		Hazard-rate	~ndvitt + watertt + size	0.768	0.271	0.060	7.411
28		Hazard-rate	~tritt	0.814	0.282	0.056	7.533
60		Hazard-rate	~ndvitt + foresttt + size	0.771	0.258	0.062	7.551
36		Hazard-rate	~observern + tritt + foresttt + size	0.643	0.269	88.048	8.266
23		Half-normal	~observern + hditt + foresttt	0.596	0.381	0.042	8.676
32		Hazard-rate	~observern + tritt + size	0.669	0.293	0.052	9.026
30		Hazard-rate	~tritt + size	0.767	0.283	0.056	9.525
50		Hazard-rate	~watertt + tritt + size	0.694	0.289	0.054	10.294
34		Hazard-rate	~tritt + foresttt + size	0.721	0.268	0.058	10.635
51		Half-normal	~watertt + foresttt + size	0.572	0.409	0.053	10.949
58		Hazard-rate	~ndvitt + tritt + size	0.649	0.269	31.127	10.992
41		Half-normal	~observern + tritt + foresttt + hditt + size	0.375	0.372	74.257	14.021
9		Half-normal	~observern	0.141	0.454	0.031	15.291
45		Half-normal	~watertt	0.130	0.455	0.031	15.848
19		Half-normal	~observern + hditt + size	0.275	0.427	0.033	16.013
1		Half-normal	~1	0.097	0.466	0.032	16.721
21		Half-normal	~observern + foresttt + size	0.146	0.446	0.031	17.019
15		Half-normal	~observern + size	0.117	0.454	0.032	17.269
47		Half-normal	~watertt + size	0.103	0.455	0.031	17.845
11		Half-normal	~foresttt	0.078	0.465	0.032	18.581
7		Half-normal	~size	0.078	0.465	0.032	18.650
37		Half-normal	~observern + tritt + foresttt	0.122	0.438	487.377	19.073
39		Half-normal	~observern + tritt	0.100	0.447	0.039	19.471
53		Half-normal	~ndvitt	0.065	0.464	0.032	20.209
5		Half-normal	~hditt	0.073	0.457	0.032	20.338
17		Half-normal	~foresttt + size	0.072	0.465	0.031	20.416
43		Half-normal	~observern + tritt + hditt + size	0.175	0.421	449.448	20.522
35		Half-normal	~observern + tritt + foresttt + size	0.135	0.438	485.884	20.829
62		Half-normal	~ndvitt + watertt + size	0.071	0.453	0.031	21.153
27		Half-normal	~tritt	0.060	0.460	0.039	21.314
31		Half-normal	~observern + tritt + size	0.085	0.447	0.039	21.398
55		Half-normal	~ndvitt + size	0.051	0.464	0.032	22.109
49		Half-normal	~watertt + tritt + size	0.080	0.448	0.039	22.139
13		Half-normal	~hditt + size	0.057	0.457	0.032	22.312
25		Half-normal	~hditt + foresttt	0.054	0.457	0.032	22.319
29		Half-normal	~tritt + size	0.052	0.460	0.039	23.188
59		Half-normal	~ndvitt + foresttt + size	0.048	0.463	0.032	23.803
33		Half-normal	~tritt + foresttt + size	0.161	0.460	537.467	24.198
57		Half-normal	~ndvitt + tritt + size	0.046	0.455	0.039	25.552

The hazard-rate in general outperforms the half-normal.

The top two models (within 2 AIC units) are the hazard-rate with both observer type, human distribution and forest as covariates, or just observer type as the sole covariate. The probability of detection within the covered area is very similar, so the differences in inference are unlikely to be meaningful.

Let’s make inferences from the top model (hazard-rate with observer type, HDI and forest as covariates).

4.3. Summary of the top model

dsmod.best <- dsmod.hr.10
  
# Output the model summary (without the abundance estimates)
print(dsmod.best$ddf)

## 
## Distance sampling analysis object
## 
## Summary for ds object
## Number of observations :  83 
## Distance range         :  0  -  725 
## AIC                    :  467.6444 
## Optimisation           :  mrds (nlminb) 
## 
## Detection function:
##  Hazard-rate key function 
## 
## Detection function parameters 
## Scale coefficient(s): 
##                   estimate        se
## (Intercept)      5.8637922 0.3526896
## observernsingle -1.1851563 0.3642999
## hditt0.02-0.03  -0.4632823 0.3834789
## hditt0.04-0.05   1.4571451 1.7099497
## hditt>0.06      -1.1302953 0.4632572
## foresttt1        1.3771032 0.6412245
## 
## Shape coefficient(s):  
##              estimate        se
## (Intercept) 0.7079553 0.2038231
## 
##                        Estimate          SE        CV
## Average p             0.2906056  0.05790504 0.1992565
## N in covered region 285.6105077 63.80722243 0.2234064

dsmod.second <- dsmod.hr.3
  
# Output the second model summary (without the abundance estimates)
print(dsmod.second$ddf)

## 
## Distance sampling analysis object
## 
## Summary for ds object
## Number of observations :  83 
## Distance range         :  0  -  725 
## AIC                    :  469.3911 
## Optimisation           :  mrds (nlminb) 
## 
## Detection function:
##  Hazard-rate key function 
## 
## Detection function parameters 
## Scale coefficient(s): 
##                   estimate        se
## (Intercept)      5.0131553 0.3578341
## observernsingle -0.6010245 0.4196622
## 
## Shape coefficient(s):  
##              estimate        se
## (Intercept) 0.4240514 0.2021686
## 
##                        Estimate          SE        CV
## Average p             0.2863336  0.05800776 0.2025880
## N in covered region 289.8716449 64.71635931 0.2232587

dsmod.obstype <- dsmod.hr.6
  
# Output the second model summary (without the abundance estimates)
print(dsmod.obstype$ddf)

## 
## Distance sampling analysis object
## 
## Summary for ds object
## Number of observations :  83 
## Distance range         :  0  -  725 
## AIC                    :  470.9706 
## Optimisation           :  mrds (nlminb) 
## 
## Detection function:
##  Hazard-rate key function 
## 
## Detection function parameters 
## Scale coefficient(s): 
##                   estimate         se
## (Intercept)      4.9668350 0.40141556
## observernsingle -0.6582242 0.41016289
## size             0.0173018 0.02785328
## 
## Shape coefficient(s):  
##              estimate       se
## (Intercept) 0.4509267 0.200856
## 
##                        Estimate          SE        CV
## Average p             0.2892537  0.05826119 0.2014190
## N in covered region 286.9453618 63.77985059 0.2222718

From the observer type model summary, we can see that the scale parameter for single observer is #negative! than that for double observer. We can also see that the effect of hdi 0.04-0.05 and forest is positive.

Next we can plot the detection function for observer type, visualising the effect of covariates also.

# Plot the overall detection function
plot(dsmod.obstype)

# Add lines for small and large herds (arbitrarily n=2 and n=20) by double observer
double.pred <- data.frame(observern=c("double", "double"), size=c(2, 15))
add_df_covar_line(dsmod.obstype, double.pred[1, ], lty=2, lwd=2, col="magenta")
add_df_covar_line(dsmod.obstype, double.pred[2, ], lty=1, lwd=2, col="magenta")

# Add the same but by single
single.pred <- data.frame(observern=c("single", "single"), size=c(2, 15))
add_df_covar_line(dsmod.obstype, single.pred[1, ], lty=2, lwd=2, col="orange")
add_df_covar_line(dsmod.obstype, single.pred[2, ], lty=1, lwd=2, col="orange")

legend("topright", legend=c("double, herd=2",
                            "double, herd=15",
                            "single, herd=2",
                            "single, herd=15"),
       title="Covariates: observer & herd size",
       lty=rep(c(2,1),times=2), lwd=2, col=rep(c("magenta", "orange"), each=2))

Now that we’ve plotted the effect of observer type and herd size, we can see that the effects are quite small in both cases. Double observer is better than single.

Let’s output a summary table of results regarding the detection function, including effective strip widths, detection probability and parameter estimates.

# Define the mean group size from the Poisson zero-truncated model we fit earlier
mean.groupsize <- round(mod.size.pred[mod.size.pred$Distance==0, "groupsize"], 0)

# Calculate the effective strip widths for double and single
esw <- predict(dsmod.obstype, data.frame(observern=c("double", "single"), size=mean.groupsize), esw=T, se=T)
esw <- data.frame(Parameter=c("ESW: double", "ESW: single"), Estimate=esw$fitted, SE=esw$se.fit)

# Extract the overall detection probability (= esw / truncation distance)
p.est <- c(summary(dsmod.obstype$ddf)$average.p, summary(dsmod.obstype$ddf)$average.p.se)
p.est <- data.frame(Parameter="Average detection prob (p)", Estimate=p.est[1], SE=p.est[2])
# Note, to get p for double and single can divide the esw by truncation distance.

# Extract the parameter estimates for the effect of vehicle and herd size
param.est <- data.frame(Parameter=c("Effect of double observer (on scale)", "Effect of herd size (on scale)"), Estimate=summary(dsmod.obstype$ddf)$coeff$key.scale[2:3, "estimate"], SE=summary(dsmod.obstype$ddf)$coeff$key.scale[2:3, "se"])

out <- rbind(esw, p.est, param.est)

kbl(out, digits=2) %>%
  kable_material(c("striped", "hover"), full_width=F) 

Parameter	Estimate	SE
ESW: double	274.07	60.90
ESW: single	160.43	45.31
Average detection prob (p)	0.29	0.06
Effect of double observer (on scale)	-0.66	0.41
Effect of herd size (on scale)	0.02	0.03

We can see that, whilst the effects of observer and herd size are substantial (as seen in the plot above), they are also highly uncertain. This is indicated by the large size of the SEs relative to the parameter estimates.Estimate: ESW for double observer is 274.07 units (meters, usually), and for single observer is 160.43 units. So, the double observer estimate is more precise.Here, 0.29 means roughly 29% of the animals present were detected, with an SE of 0.06. Negative -0.66 suggests that the “scale” parameter decreases with double observers, but the SE is 0.41, so the effect is not highly precise (likely not significant).Estimate 0.02 is very small with SE 0.03, meaning herd size has almost no effect on detection probability in this dataset.

5. Excluding rugged areas from the survey area

Now that we have a best model, we can make inferences about density. However, if we want to estimate abundance in the survey area, we need to first exclude the rugged mountain areas which are likely non-habitat for gazelles or marginal habitat at best.

We can exclude the rugged areas by finding a threshold value for the terrain ruggedness index which more-or-less defines the marginal habitat for gazelles.

There are model-based approaches we could use to guide this (e.g. using a Poisson generalized linear model of the counts along the transects), but we would still have to define the threshold using an arbitrary cutoff.

Here we will simply plot the range of ruggedness values in the study area and visualise what range of values are: i) covered by the transects, and ii) covered by the gazelle detections.

# Get TRI values within the survey area
tri.vals <- values(mask(tri, st_as_sf(survey.area)))

# Get TRI along the transects (every 200 m)
transect.points <- st_as_sf(st_cast(st_zm(st_line_sample (transects, density=1/100), drop=T), "POINT"))
tri.transects <- extract(tri, transect.points)

# Get TRI for the detections
tri.detections <- extract(tri, object.locations[object.locations$species=="goitered_gazelle", ])

ggplot(data=data.frame(tri.vals), aes(x=tri.vals)) +
  geom_histogram(aes(y=after_stat(density)), fill="grey", alpha=0.7) +
  geom_rug(data=data.frame(tri.transects), aes(x=tri.transects), sides="t") +
  geom_rug(data=data.frame(tri.detections), aes(x=tri.detections), sides="b", colour="red") +
  labs(x="Terrain ruggedness index") +
  theme_minimal()

The transects sometimes were mostly in areas of TRI < 0.05, but sometimes went into more rugged terrain. 11 gazelle detections were made in these areas, suggesting it is marginal habitat (>0.05).

Let’s threshold the ruggedness map with a value of 0.05 and visually inspect the result.

tri.thresh <- 0.05

# Threshold the tri raster
mountains <- app(rast(tri), function(x) {
  x[x < tri.thresh] <- NA
  x[x >= tri.thresh] <- 1
  x
})

# Convert the mountain raster to a polygon (this looks a bit better than the raster when plotted)
mountains.sf <- st_as_sf(as.polygons(mountains))

# Plot it
leaflet() %>% 
  addTiles(urlTemplate = "http://mt1.google.com/vt/lyrs=s&x={x}&y={y}&z={z}") %>% 
  #addRasterImage(mountains, color="white", opacity=0.95) %>%
  addPolygons(data=st_transform(mountains.sf, "EPSG:4326"), color="white", weight=1.3, opacity=0.35, fillColor="white", fillOpacity=0.7) %>%
  addPolylines(data=st_transform(survey.area, "EPSG:4326"), color="white", weight=1.3, opacity=0.6, dashArray="4") %>%
  addPolylines(data=st_transform(st_zm(transects.simple, drop=T, what="ZM"), "EPSG:4326"), color="black", weight=1.5, opacity=1) %>%
  addCircleMarkers(data=st_transform(object.locations[object.locations$species=="goitered_gazelle", ], "EPSG:4326"), radius=3, stroke=F, fillOpacity=0.8, color="orange")

Using a threshold of 0.05 actually does not exclude all of the mountain areas, only the highest and most rugged areas. Given that the gazelles are sometimes found in the mountain areas, especially in the lower slopes, this might actually be a realistic picture of gazelle habitat. We will proceed on this basis for now.

6. Estimating density and abundance

We are now ready to estimate density and abundance across the three surveys.

# We first adjust the areas in the region.table to exclude the mountain areas

## Calculate areas covered by the rugged terrain
tost.mtn.area <- as.numeric(st_area(st_intersection(mountains.sf, tost.area)) / 1000^2)
#ggn.mtn.area <- as.numeric(st_area(st_intersection(mountains.sf, ggn.area)) / 1000^2)

## Subtract them from the overall area
region.table.updated <- region.table
region.table.updated$Area <- c(region.table.updated$Area[1] - tost.mtn.area ) 

# And then update the model with these new areas
dsmod.best.updated <- ds(obs.table, formula=~observern+hditt+foresttt, key="hr", adjustment=NULL, cutpoints=ds.cutpoints, monotonicity="none", region_table=region.table.updated, sample_table=sample.table, obs_table=obs.table, convert_units=0.001)

# Output the model summary
dsmod.best.updated$dht

## Abundance and density estimates from distance sampling
## Variance       : R2, N/L 
## 
## Summary statistics
## 
##      Region     Area CoveredArea   Effort  n  k         ER       se.ER
## 1 Tost-2025 15987.13    1921.561 1325.215 83 28 0.06263137 0.008534237
##       cv.ER
## 1 0.1362614
## 
## Summary for clusters
## 
## Abundance:
##   Region Estimate      se        cv      lcl      ucl       df
## 1  Total  2376.24 644.654 0.2712916 1397.795 4039.589 77.22753
## 
## Density:
##   Region  Estimate         se        cv        lcl       ucl       df
## 1  Total 0.1486346 0.04032332 0.2712916 0.08743255 0.2526777 77.22753
## 
## Summary for individuals
## 
## Abundance:
##   Region Estimate       se        cv     lcl      ucl     df
## 1  Total 14287.01 4024.049 0.2816579 8226.18 24813.31 63.659
## 
## Density:
##   Region  Estimate        se        cv       lcl      ucl     df
## 1  Total 0.8936574 0.2517056 0.2816579 0.5145503 1.552081 63.659
## 
## Expected cluster size
##   Region Expected.S se.Expected.S cv.Expected.S
## 1  Total   6.012444     0.7450291     0.1239145

The coefficients of variation are rather large, indicating that there is a lot of uncertainty in the estimates. An ideal target for monitoring is a CV < 0.2. CVs of > 0.20 are particularly okey.

Finally, we can plot the density and abundance estimates.

# Extract the density and abundance summary tables
density <- data.frame(Variable="Density (per km2)", summary(dsmod.best.updated)$dht$individuals$D)
abundance <- data.frame(Variable="Abundance", summary(dsmod.best.updated)$dht$individuals$N)

# Combine them
final.estimates <- rbind(density, abundance)

# Put density first
final.estimates$Variable <- factor(final.estimates$Variable, levels=c("Density (per km2)", "Abundance"))

# Remove the total abundance, as it doesn't make sense
final.estimates <- final.estimates[!(final.estimates$Label=="Total" & final.estimates$Variable=="Average"), ]

# Change Total to "Average"
final.estimates$Label[final.estimates$Label=="Total"] <- "Average"

# Get a silhouette of a gazelle to make plot more interesting
gazelle.img <- get_phylopic(uuid=get_uuid(name = "Gazella gazella", n = 2))

# Simple plot without gazella silhouette  
#ggplot(final.estimates, aes(x=Label, y=Estimate, colour=Label)) +
 # geom_point(size=2) +
  #geom_errorbar(aes(ymin=lcl, ymax=ucl), width=0.2) +
  #theme_minimal() +
  #theme(legend.position="none", axis.title=element_blank()) +
  #facet_wrap(~Variable, scales="free")

# Plot the density estimates
plot.dens <- ggplot(final.estimates[final.estimates$Variable=="Density (per km2)", ], aes(x=Label, y=Estimate, colour=Label)) +
  geom_point(size=2) +
  geom_errorbar(aes(ymin=lcl, ymax=ucl), width=0.1) +
 labs(y=expression(paste("Density (animals per km"^"2", ")"))) +
 lims(y=c(0, 1.6)) +
   add_phylopic(img=gazelle.img, x=3.2, y=0.3, height=0.8, alpha=0.3) +
  theme_minimal() +
  theme(legend.position="none", axis.title.x=element_blank()) 

# Plot the abundance estimates
plot.abund <- ggplot(final.estimates[final.estimates$Variable=="Abundance", ], aes(x=Label, y=Estimate, colour=Label)) +
  geom_point(size=2) +
  geom_errorbar(aes(ymin=lcl, ymax=ucl), width=0.1) +
  labs(y="Estimated abundance") +
  lims(y=c(0, 25000)) +
 add_phylopic(img=gazelle.img, x=2.0, y=3000, height=11000, alpha=0.3) +
  theme_minimal() +
  theme(legend.position="none", axis.title.x=element_blank()) 

plot.dens + plot.abund

# Spatial Density (Нягтрал)

library(ggplot2)
library(dplyr)

# 1. Координат багануудыг шалгах
if(!all(c("objectX", "objectY") %in% names(dat.ggaz))) {
  stop("objectX болон objectY баганууд датанд байхгүй байна!")
}

# 2. NA-гүй координат бэлдэх
coords <- dat.ggaz %>%
  select(objectX, objectY) %>%
  na.omit()

# 3. 2D KDE зураг
p_spatial <- ggplot(coords, aes(x = objectX, y = objectY)) +
  stat_density_2d(aes(fill = ..level..), geom = "polygon", color = "white") +
  scale_fill_viridis_c(option = "magma") +
  coord_equal() +
  labs(
    title = "Spatial Kernel Density Map (2D KDE)",
    x = "objectX",
    y = "objectY",
    fill = "Density"
  ) +
  theme_minimal(base_size = 14)

# RMarkdown-д гаргах
print(p_spatial)

# Хадгалах (заавал биш)
ggsave("outputs/spatial_KDE_map.png", p_spatial,
       width = 10, height = 6, dpi = 200)

7. Conclusions

• This is a different result from the combined analysis of data from the two studies, the Tost and the GGN. Otherwise, the analysis model and the pattern of the effects of these covariates appear to be different in the two study areas. This analysis shows that the confidence variation has a high value (>0.20). In general, observer type appears to be a key influencing variable, although the statistical precision is unclear. The number of detections in the study is small, with perfect detection at 0-50 m, a sharp drop at 100-400 m, and a severe lack of detection at 400-2000 m, which significantly affects the detection function. This may be due to open ground visibility, animal behavior, observer error, and binocular delay.

• There is a strong and statistically clear effect of observer type ( unclear effect HDI and forest) on detection, with double observer(with 2 rangers)–based surveys having a larger detection width

• Overall density across the this surveys was 0.89 individuals per km²

• By buffering the transects by 8 km, we can define a large survey area, in which we estimate there are 14,287 goitered gazelles in Gobi gurvansaikhan NP (around 16000 km²).

###Rdatafile
#save(data,file="Tost_Gazelle_DS_2025_v4_dataframe.Rda")
#load(file="Tost_Gazelle_DS_2025_v4_dataframe.Rda")
#save.image("Tost_Gazelle_DS_2025_v4")
#save.image("Tost_Gazelle_DS_2025_v4_workspace.Rdata")